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Abstract 

Linear viscoelastic properties for a dilute polymer solution are predicted by modeling the solution 
as a suspension of non-interacting bead-spring chains. The present model, unlike the Rouse model, can 
describe the solution's rheological behavior even when the solvent quality is good, since excluded volume 
effects are explicitly taken into account through a narrow Gaussian repulsive potential between pairs 
of beads in a bead-spring chain. The use of the narrow Gaussian potential, which tends to the more 
commonly used 5-function repulsive potential in the limit of a width parameter d going to zero, enables 
the performance of Brownian dynamics simulations. The simulations results, which describe the exact 
behavior of the model, indicate that for chains of arbitrary but finite length, a 5-function potential leads 
to equilibrium and zero shear rate properties which are identical to the predictions of the Rouse model. 
On the other hand, a non-zero value of d gives rise to a prediction of swelling at equilibrium, and an 
increase in zero shear rate properties relative to their Rouse model values. The use of a 5-function 
potential appears to be justified in the limit of infinite chain length. The exact simulation results are 
compared with those obtained with an approximate solution which is based on the assumption that the 
non-equilibrium configurational distribution function is Gaussian. The Gaussian approximation is shown 
to be exact to first order in the strength of excluded volume interaction, and is found to be accurate 
above a threshold value of d, for given values of chain length and strength of excluded volume interaction. 



*Present address: Department of Chemical Engineering, Monash University, Victoria 3800, Australia 



1 Introduction 



The simplest model, within the context of Polymer Kinetic Theory, to describe the rheological behavior of 
dilute polymer solutions is the Rouse model.El The Rouse model represents the macromolecule by a linear 
chain of identical beads connected by Hookean springs, and assumes that the solvent influences the motion 
of the beads by exerting a drag force and a Brownian force. While the Rouse model is able to explain the 
existence of viscoelasticity in polymer solutions by predicting a constant non-zero first normal stress difference 
in simple shear flow, it cannot predict several other features of dilute solution behavior, such as the existence 
of a nonzero second normal stress difference, the existence of shear rate dependent viscometric functions, or 
the correct molecular weight dependence of material functions. Over the past decade, considerable progress, 
has been made by incorporating the effect of fluctuating hydrodynamic interaction into the Rouse model. 
These models are able to predict the molecular weight dependence of the material functions accurately. 
They also predict a nonzero second normal stress difference, and shear rate dependent viscometric functions. 
However, since they neglect the existence of excluded volume interactions among parts of the polymer chain, 
they are strictly applicable only to theta solutions. 

Recently, Prakash and OttingerQ examined the influence of excluded volume effects on the rheological 
behavior of dilute polymer solutions by representing the polymer molecule with a Hookean dumbbell model, 
and using a narrow Gaussian repulsive potential to describe the excluded volume interactions between the 
beads of the dumbbell. The narrow Gaussian potential tends to a (^-function potential in the limit of a 
parameter, that describes the width of the potential, going to zero. It can therefore be used to evaluate 
results obtained with the singular S- function potential. It was shown by them that the use of rf-i-4- function 
potential between the beads, which is commonly used in static theories for polymer solutions,Er£3 leads to 
no change in the equilibrium or dynamic properties of the dumbbell when compared to the case where no 
excluded volume interactions are taken into account. They also found that assuming that the non-equilibrium 
configurational distribution function is a Gaussian leads to accurate predictions of viscometric functions in 
a certain range of parameter values. These results suggest that it would be worthwhile examining longer 
bead-spring chains. Firstly, it is interesting to see if the problem with the <5-function potential can be resolved 
when there are more beads in the bead-spring chain. Secondly, it is important to find out if the Gaussian 
approximation is accurate even for longer chains. The purpose of this paper is to attempt to answer these 
questions, in the linear viscoelastic limit, by extending the methodology developed in the earlier paper to 
the case of bead-spring chains. The same issues, in the context of steady shear flows at finite shear rates, 
will be addressed in a subsequent paper. 

As in the dumbbell paper, we confine our attention to excluded volume interactions alone, and neglect the 
presence of hydrodynamic interactions. This clearly implies — since it is essential to include hydrodynamic 
interaction effects for a proper description of the dynamic behavior of dilute solutions — that the results of the 
present paper are not yet directly comparable with experiments. They represent a preliminary step in that 
direction. It is felt that the inclusion of hydrodynamic interaction would make the theory significantly more 
complex before the role of excluded volume interactions is properly understood. The aim of this work is to 
develop and carefully evaluate the Gaussian approximation for excluded volume interactions. The Gaussian 
approximation has already been shown to be excellent for the treatment of hydrodynamic interaction effects. □ 
If it also proves to be accurate for the treatment of excluded volume effects, then it would be extremely useful 
for describing the combined effects of hydrodynamic interaction and excluded volume. It should be noted 
that, in principle, the development of such an approximation does not pose any fundamental problems. 

This paper is organized as follows. In the next section, the basic equations governing the dynamics of 
Rouse chains with excluded volume interactions are discussed. A retarded motion expansion for the stress 
tensor is derived in section 3, and exact expressions for the zero shear rate viscometric functions in simple 
shear flow are obtained. The implications of these results for a ^-function excluded volume potential are 
then discussed. In section 4, the Brownian dynamics simulation algorithm used in this work is described. 
Section 5 is devoted to the development of the Gaussian approximation for the configurational distribution 
function. Exact expressions for linear viscoelastic properties are derived through a codeformational memory- 
integral expansion. In section 6, a first order perturbation expansion in the strength of the excluded volume 
interaction is carried out. This proves to be very helpful in understanding the nature of the Gaussian 
approximation. The results of the various exact and approximate treatments are compared and discussed in 
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section 7, and the main conclusions of the paper are summarized in section 8. 



2 Basic Equations 

The instantaneous configuration of a linear bead-spring chain, which consists of N beads connected together 
by (N — 1) Hookean springs, is specified by the bead position vectors r u , v = 1, 2, . . . , N, in a laboratory- 
fixed coordinate system. The Newtonian solvent, in which the chain is suspended, is assumed to have a 
homogeneous velocity field — that is, at position r and time t, the velocity is given by v(r, t) = Vo + n(t) ■ r, 
where vq is a constant vector and n(t) is a traceless tensor. 

The microscopic picture of the intra-molecular forces within the bead-spring chain is one in which the 
presence of excluded volume interactions between the beads causes the chain to swell, while on the other 
hand, the entropic retractive force of the springs draws the beads together and opposes the excluded volume 
force. This is modeled by writing the potential energy (f> of the bead-spring chain as a sum of the potential 
energy of the springs S 1 , and the potential energy due to excluded volume interactions E. The potential 
energy S is the sum of the potential energies of all the springs in the chain, and is given by, 

1 JV_1 

s = -Hj2Qi-Qi (!) 

i=l 

where, H is the spring constant, and Q i = r%+i — T{, is the bead connector vector between the beads i and 
i + 1. The excluded volume potential energy E is found by summing the interaction energy over all pairs of 
beads \i and v, 

1 N 

E= g J2 E(r v -r M ) (2) 

where, E (r„ — r^) is a short-range function. It is usually assumed to be a <5-function potential in static 
theories for polymer solutions, 

E(r v -r lt )=vk B T6(r„-r li ) (3) 

where, v is the excluded volume parameter (with dimensions of volume), /cb is Boltzmann's constant, and T is 
the absolute temperature. In this work, we regularise the S- function potential, and assume that E [r v — tv) 
is given by a narrow Gaussian potential, 

^->(-W) 

where, d is a parameter that describes the width of the potential (it represents, in some sense, the extent 
of excluded volume interactions), and r„ M = r„ — r M , is the vector between beads /x and v. In the limit d 
tending to zero, the narrow Gaussian potential becomes a ^-function potential. 

The intra-molecular force on a bead is, — —(dtp/dr^), can consequently be written as, F)f^ = 



Fl s) +F[ e \ where, 



OTv k=i 



dE * d 



dr v ^-f dr v 



In eq ^, is an (N -l)xJV matrix defined by, Bk v = Sk+i, v — $kv, with 5k u denoting the Kronecker 
delta. 

For homogeneous flows, the internal configurations of the bead-spring chain are expected to be indepen- 
dent of the location of the centre of mass. Consequently, it is assumed that the configurational distribution 
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function ip depends only on the (N — 1) bead connector vectors Q k . The diffusion equation that governs 
ip (Qi, ■ ■ ■ , Q N _i, t), for a system with an intra-molecular potential energy <fr as described above, can then 
be shown to be given by, 

rj , N-l a / tt N ~ 1 I N \ 

where, £ is the bead friction coefficient (which, for spherical beads with radius a, in a solvent with viscosity 
r] s , is given by the Stokes expression: £ = 6irr] s a), and Ajk is the Rouse matrix, 

n ( 2 for |j - fc| = 0, 

^ife = £ BjvBkv = I -1 for \j -k\ = l, (8) 
" — v otherwise 

The time evolution of the average of any arbitrary quantity, carried out with the configurational distri- 
bution function if), can be obtained from the diffusion equation. In particular, by multiplying eq|7|by QjQ^, 
and integrating over all configurations, the following time evolution equation for the second moments of the 
bead connector vectors is obtained, 

j t {QjQi) = K ■ (QjQk) + (QjQk) ■ kT + A jk l 

H 

- 7 E { (Q 3 Qm) Arnk + A jm (Q m Q k ) } + Yjk (9) 
' m=l 

where, 1 is the unit tensor, and, 

1 N 

^= E {(Q 3 F^)B Hl+ B m (F^Q k )} (10) 

The term Yjk, which arises due to the presence of excluded volume interactions, does not appear in the 
second moment equation for the Rouse model. Due to this term, which in general involves higher order 
moments, eq ^, is not a closed equation for (QjQ k ). As will be discussed in greater detail in the section on 
the Gaussian approximation, finding an approximate solution for the present model involves making eq [)] a 
closed equation for the second moments. 

The polymer contribution to the stress tensor — for models with-arbitrary intra-molecular potential forces 
but no internal constraints — is given by the Kramers expression,^ 

N-l 

r'P = -n p H (QkQk) + Z + (N-1) n p k B T 1 (11) 



fc=i 



where, 



N N-l 



Z = « p ^^^(Q t ^) (12) 



v=\ fe=l 



Here, n p is the number density of polymers, and B vk is a JVx (N— 1) matrix defined by, B vk = k/N—Q (k—v), 
with (A; — v) denoting a Heaviside step function. 

It is clear from cq |ll| that there are two reasons why the presence of excluded volume interactions leads 
to a stress tensor that is different from that obtained in the Rouse model. Firstly, there is an additional 
term represented by Z which is the direct influence of excluded volume effects. Secondly, there is an indirect 
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influence due to a change in the contribution of the term ^ifclTi 1 (QkQk)j relative to its contribution in the 
Rouse case. For a (5-funcijon excluded volume potential, it can be shown that the direct contribution to the 
stress tensor is isotropic.S On the other hand, for the narrow Gaussian potential, Z is not isotropic unless d 
is equal to zero. It is therefore important to use the complete form of the Kramers expression, eq [ll], when 
carrying out simulations with an excluded volume potential that is not a (5-function potential. 

All the rheological properties of interest can be obtained once the stress tensor in eq [ll] is evaluated. In 
the next section, a retarded motion expansion for the stress tensor is derived. 



3 Retarded Motion Expansion 

A retarded motion expansion for the stress tensor can be obtained by extending the derivation carried out 
previously for the dumbbell modefl to the case of bead-spring chains. The dumbbell model derivation was, 
in turn, an adaptation of a similar development for the FENE dumbbell model.EU The argument in all these 
cases rests basically on seeking a solution of the diffusion equation, eq 0, of the following form, 

i>(Ql, Qn-V*) = V>eq(Ql, • ■ ■ , Qn-i) 0fl (Qlj • • • , Qn-V*) ( 13 ) 

where, tp cq is the equilibrium distribution function given by, 

A q (Qi,---,Q N -i)=Mc q e~^ kBT (14) 

with A/" oq denoting the normalization constant, and 0a is the correction to ip cq due to flow — appropriately 
termed the flow contribution. 

The governing partial differential equation for 0a (Q l7 . . . , Q N _ ± , t) can be obtained by substituting eq|l3| 
into the diffusion equation, eq ^. It turns out that, regardless of the form of the excluded volume potential, 
at steady state, an exact solution to this partial differential equation can be found for all homogeneous 
potential flows. For more general homogeneous flows, however, one can only obtain a perturbative solution 
of the form, 

fl (Q l5 . . . , Q N -!,t) = 1 + 01 + 02 + 03 + ■ ■ • (15) 

where 0fc is of order fc in the velocity gradient. 

Partial differential equations governing each of the 0^ may be derived by substituting eq [l5] into the 
partial differential equation for 0a and equating terms of like order. The forms of the functions 0& can then 
be guessed by requiring that they fulfill certain conditions.til In the present instance, we only find the form 
of 0i, since our interest is confined to zero shear rate properties. One can show that, 



jV-1 

C m nQ m -i-Qn (16) 



ik B T ^ 

where, 7 is the rate of strain tensor, 7 = Vi> + Vv T , and C mn is the Kramers matrix. The Kramers matrix 
is the inverse of the Rouse matrix, and is defined by, 

JY 

C mn = B urn B vn = min (m, n) - mn/N (17) 
i^=i 

In order to proceed further, we need to show that the present model satisfies the Giesekus expression for 
the stress tensor Jiil Upon multiplying eq^ with Y^ m ~n=i Cmn Q m Qm an d integrating over all configurations, 
we can show that, 

d N 1 N 1 

m,n— 1 m,n— 1 

91. T Off N ~ 1 N Ar " 1 

= ^(«-l)l-yE (Q m Qm) + EE B um (Q m FW) (18) 
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On combining this equation with eq[Tl| for the stress tensor, it is straight forward to see that the Giesekus 
expression is indeed satisfied. At steady state the Giesekus expression reduces to, 



2 



N-l 



]T C mn {K-(Q m Q n ) + (Q m Q n )-K T } 



(19) 



Clearly, the stress tensor at steady state can be found once the average (Q m Q n ) is evaluated. This can 
be done, correct to first order in velocity gradients, by using the power series expansion for ipa, eq 15, with 
the specific form for <f>\ in eq [l6| The following retarded motion expansion for the stress tensor, correct to 
second order in velocity gradients and valid for arbitrary homogeneous flows, is then obtained, 
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JV-l 



« • (Q m Qn) cq + (QmQn) ca ■ « 



n P C 2 



N-l N-l 



K ■ (Q m Q n (Qj ' 7 • Qk)) e 



m,n— 1 j,k— 1 

;(Q,-7-Q fc )Q m Q„) ca -« T 



(20) 



where, denotes the average of any arbitrary quantity X with the equilibrium distribution function 

One can see clearly from eq[2C] that rheological properties, at small values of the velocity gradient, can 
be obtained by merely evaluating equilibrium averages. The special case of steady simple shear flow in the 
limit of zero shear rate is considered below. 



3.1 Zero Shear Rate Viscometric Functions 

Steady simple shear flows are described by a tensor n which has the following matrix representation in the 
laboratory-fixed coordinate system, 

/0 1 0\ 

k = 7 (21) 
\0 0/ 

where 7 is the constant shear rate. 

The three independent material functions used to characterize such flows are the viscosity, rj p , and the 
first and second normal stress difference coefficients, and ^2, respectively. These functions are defined 
by the following relations ,E3 

r^ = -iv P ; t-L-t£, = - 7 2 *i; r^-rf z = -7 2 * 2 (22) 

The components of the stress tensor in simple shear flow, for small values of the shear rate 7, can be 
found by substituting eq [2l] for the rate of strain tensor, into eq ^0|. This leads to, 

t . N-l r2 ■ 2 N ~ 1 N ~ 1 

p _^p_47 v (y y \ - rip ^ 7 V Vr r (y y x y\ 

l xy — 2 rnn \ m ™/cq 4fc B T Zv °mn L 'M\- I ni J n A P I 9/ cq 

m,n— 1 m,n— 1 p : q—l 



N-i ^ C 2 7 2 N ~ X N ~ 1 

j.i ~ — "'pC'7 ^ ] Cran (■^m^i) eq 2k^T ^ < ^ t C mn Cp q (X m Y n X p Y q ^ ^ 



"p C 7 
2fc B T 

m,n— 1 rn,n— 1 p,q=l 

= ^ = (23) 
where, (X m ,17„,Z m ) are the Cartesian components of the bead connector vector Q m . 
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Using the symmetry property of the potential energy <f>, which remains unchanged when the sign of the 
Yk component of all the bead connector vectors Q k ; k — 1,2, ... , (N — 1), is reversed, we can show that, 
(X m Y n } e — 0. From the definitions of the viscometric functions in eq ^2], it is straight forward to show that, 
in the limit of zero shear rate, the following exact expressions for the zero shear rate viscometric functions 
are obtained. 



c N-l 

^ E^niQm-QnU (24) 



*i,o = ^ E E C ^ C vi (X m Y n X p Y q ) cq (25) 

m.n—1 p,q=X 

* 2 ,o = (26) 

In order to derive eq |24|, we have used the fact that, since 4> is the same function of X p , Y p , and Z p for all p, 
(X p X q ) cq = (YpYq} — (Z p Z q ) cq . Equation p6| indicates that the inclusion of excluded volume interactions 
alone is not sufficient to lead to the prediction of a non-zero second normal stress difference. The proper 
inclusion of hydrodynamic interaction is required. 

It is interesting ta note, by making use of eq |24|, that the mean square radius of gyration at equilibrium, 
which is defined as, til 

1 N f 

( R l > cq = n £ d Qi d Q2 ■ ■ ■ dQjv-i (r* - r c ) ■ (r„ - r c ) ^ (27) 

u=l J 

(where, r c is the position of the center of mass), is related to the zero shear rate viscosity by, 

V P ,o = ^N(R* g ) eq (28) 

An alternative expression for the zero shear rate viscosity, which will prove very useful subsequently, can 
also be obtained from eq p3, 

r N 

i 



N-l N-l 



In order to derive eqs 28 and |2S|, equations which relate the bead connector vector coordinates to bead 
position vector coordinates, summarized for example, in Chapter 11 and Table 15.1-1 of Chapter 15 of the 
text book by Bird et al.,EJ have been used. 

The evaluation of the equilibrium averages in eqs ^5] and for various values of the parameters in the 
narrow Gaussian potential, and for various chain lengths N, have been carried out here with the help of 
Brownian dynamics simulations. More details of these simulations are given subsequently. In the special 
case of the extent of excluded volume interactions d going to zero or infinity, we had shown earlier for a 
Hookean dumbbell model that the values of r/p.or^nd ^i,o remain unchanged from the values that they have 
in the absence of excluded volume interactions .0 In the next section, we consider the same limits for the 
more general case of bead-spring chains of arbitrary (but finite) length. 

3.2 The Limits d — > and d — > oo 

The average in eq ^ can be evaluated with the distribution function ■0 eq (Q 1 , . . . , Q^_ 1 ), or equivalently, 
with the distribution function Peq^w/i)) which is a contracted distribution function for each vector r vpil and 
which is defined by, 



P e q(r vll ) = JdQ 1 dQ 2 . . . dQ N _ x S - E Q 3 ■ j V'cq 



(30) 



We have assumed here, without loss of generality, that v > 
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In the Rouse model, as is well known, the equilibrium distribution function is Gaussian, 

<(Ox,...,o J ,_ 1 )= n (db) 7 cxp (-2^t^-^) (3l) 

A superscript or subscript 'i?' on any quantity will henceforth indicate a quantity defined or evaluated in 
the Rouse model. The distribution function P^fc-^) can then be evaluated, by substituting eq|3l] and the 
Fourier representation of a (5-function, into eq pOUlj 

< <r "-> " ( frtrVd ) ' oxp (-W^ < 32 > 

The absolute value \v — fi\ indicates that this expression is valid regardless of whether v or fi is greater. This 
is another well known result of the Rouse model, namely, at equilibrium, the vector r„ M between any two 
beads y, and ^, also obeys a Gaussian distribution. 

A similar procedure can be adopted to evaluate P eq (r 1/At ), in the presence of excluded volume interactions, 
by substituting eq [l4| and the Fourier representation of a (5-function, into eq 30. We show in appendix |A| 
that, 

lim P cq (rv) = P* (tv) (33) 



d-*0 

, d^c 



As a result, for all quantities X(r„ M ), such that the product X(r vll ) P eq ( , r l//i ) remains bounded for all r„ M , 
lim g_, (X(r vil )} = (X(r Vft ))\ It follows from eq|| that, 



lim ry p , = r)f Q (34) 

d^O 



Thus, the polymer contribution to the viscosities in the limit of zero shear rate, for chains of arbitrary (but 
finite) length, in (i) the presence of (5-function excluded volume interactions, and (ii) the absence of excluded 
volume interactions (the Rouse model), are identical to each other. Brownian dynamics simulations, details 
of which are given in the section below, indicate that this is also true for the first normal stress difference 
coefficients. 



4 Brownian Dynamics Simulations 

The equilibrium averages in eqs [2^ and |2^, as mentioned above, can be evaluated with the help of Brownian 
dynamics simulations. As a result, exact numerical predictions of the zero shear rate viscometric functions 
can be obtained. Brownian dynamics simulations basically involve the numerical solution of the-Jto stochastic 
differential equation that corresponds to the diffusion equation, eq[7| Using standard methodaia to transcribe 
a Fokker-Planck equation to a stochastic differential equation, one can show that eq [?] is equivalent to the 
following system of (N — 1) stochastic differential equations for the connector vectors Qj, 

dQ 3 =^ K . Qj -^A jk ^dt + j2 dW v (35) 

where, W u is a 3N dimensional Wiener process. ,_. 

A second order predictor-corrector algorithm with time-step extrapolatiorJ13 was used for the numerical 
solution of eq [5f| Steady-state expectations at equilibrium were obtained by setting n = 0, and simulating 
a single long trajectory. This is justified based on the assumption of ergodicity.Ej 
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5 The Gaussian Approximation 



A crucial step in the calculation of the rheological properties predicted by the present model is the evalua- 
tion of the complex moments that occur in Kramers expression. The Gaussian approximation — which has 
prevkxusbz-been shown to be useful in the treatment of hydrodynamic interaction and internal viscosity ef- 
fectsBli-jEj — consists essentially of reducing complex higher order moments to functions of only second order 
moments by assuming that the non-equilibrium configurational distribution function is a Gaussian distribu- 
tion, and subsequently, evaluating these second order moments by integrating a time evolution equation. 

For the narrow Gaussian potential, the complex moment (Q fe f 1 ^ B ' ) ), which appears in the quantity Z 
on the right hand side of Kramers expression, eq [H], can be rewritten in terms of averages of the form: 
(Q k Q n E (r v — r M ) \. Assuming that tp is a Gaussian distribution of the form, 

V (Q x , ■ ■ • , Q N -i,t) = Af(t) exp [ - ~ Q ■ (Oik ' Q k ] (36) 

j, k 

where, the (N — 1) x (A — 1) matrix of tensor components tTjk (with crjk = (QjQk) an d (Tjk — ojy) uniquely 
characterizes the Gaussian distribution and M(t) is.the normalization factor, and using general decomposition 
rules for the moments of a Gaussian distribution,!] one can show that, 

(Q m Q n E(r u -r,)) r/, ' uT ' 



( 2 ^ 3/2 ^det ([J 2 1 + <r„ M r VM )]^ 

{(Q ro Q») - (Q m r vv ) ■ [d 2 1 + (r Vfl r Ufl )] 1 • {r vll Q n )\ (37) 

The vector r ufi also obeys a Gaussian distribution since it is a sum of Gaussian distributed bead connector 
vectors. As a result, the right hand side of eq ^ involves only second moments, and averages which can be 
evaluated by Gaussian integrals. 

In the Gaussian approximation therefore, Kramers expression for the stress tensor can be rewritten as, 

N-l 

T p = -n p H^2a kk + Z+(N~l)n p k B Tl (38) 
fc=i 

where, 



1 N 

Z = -z*n p k B T ^ o^M" n (^M) ( 39 ) 



V, ,1=1 



In eq 39, the function II (o^) is given by 



n(a-^) = 1 (40) 

/det ([d* 2 1 + dup]) 



with the tensors b vfl defined by 



; T 



max(/i,i^)-l 

ti 



k B T 

j.k — mm(p..v) 



The quantities z* and d* are non-dimensional versions of the two parameters, v and d, which characterize 
the narrow Gaussian potential. They are defined by, 



'-'Ve? (42) 
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While z* measures the strength of the excluded volume interaction, d* is a measure of the extent of excluded 
volume interaction. 

In the limit of d* — > 0, it is straight forward to see that the tensor Z becomes isotropic. As a result, the 
direct contribution to the stress tensor has no influence on the rheological properties of the polymer solution 
only when a 5-function potential is used to represent excluded volume interactions. 

All that remains to be done in order to evaluate the stress tensor is to find the components of the 
covariance matrix (Tjk- A system of 9 (N — l) 2 coupled ordinary differential equations for oj/- can be obtained 
from the time evolution equation for the second moments, eq ^. As mentioned earlier, in the presence of 
excluded volume interactions, eq |^ also involves higher order moments due to the occurrence of the term 
Yjk, and consequently, it is not in general a closed equation for the second moments. However, these higher 
order moments can also be reduced to second order moments with the help of the decomposition result, 
eq |3^. In the Gaussian approximation, the second moment equation can therefore be rewritten as, 

(Tjk = K ■ (Tjk + (Tjk ■ K T H ~ Ajk 1 — ^ [(Tj m A m k + Aj m CTmk] + Y jk (43) 



where, 



Y jk = z* ( — j ^ Wjm ■ Afem + A jm ■ CT mk ] (44) 



In eq|44|, the (N — 1) x (N — 1) matrix of tensor components A jm is defined by, 



N 



^■3 m = (Bj+i, m — B^ m ) n(o- J+ i iAl ) — (Bj m — B^ m ) n(oj^) > (45) 
11=1 ^ J 

For any homogeneous flow, rheological properties predicted by the Gaussian approximation can be ob- 
tained by appropriately choosing the tensor k, solving the differential equations, eqs |]| for Ojk, and sub- 
stituting the result into Kramers expression, eq Ej8[ In this paper, we confine attention to the prediction of 
linear viscoelastic properties, namely material functions in small amplitude oscillatory shear flow, and zero 
shear rate viscometric functions. 

Linear viscoelastic properties predicted by the Gaussian approximation can be obtained by constructing 
a codeformational memory-integral expansion. This is done by expanding the tensors (Tjk, in terms of 
deviations from their isotropic equilibrium solution, up to first order in velocity gradient, 

(Tjk = fjk 1 + e jk + ... (46) 

where, the tensors fjk 1 represent equilibrium second moments in the Gaussian approximation, and the 
tensors ejk are the first order corrections. Since the details of the calculation are not very illuminating, they 
are given in appendix |b| and only the results are summarized below. 

The first order codeformational memory-integral expansion derived by the above procedure has the form, 

dsG(t-s) 1[1] (t, S ) (47) 

where, 7nj is the codeformational rate-of-strain tensor,El and the memory function G(t) is given by eq |85| 
in appendix [b]. This expression can now be used to obtain exact expressions for material functions in small 
amplitude oscillatory shear flow, and for the zero shear rate viscosity and first normal stress difference 
coefficient in steady shear flow, as shown below. 

Small amplitude oscillatory shear flow is characterized by a tensor n(t) given by 

/0 1 0\ 

K {t) = 70 cos wt (48) 
I / 
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where, 70 is the amplitude, and to is the frequency of oscillationsJn the plane of flow. The yx component of 
the polymer contribution to the shear stress is then defined by ,£3 



T yx — — 7o cos — 7]"(uj) 70 sin tot (49) 

where 77' and 77" are the material functions characterizing oscillatory shear flow. They can be represented in 
a combined form as the complex viscosity, 77* = rf — if]", and 77* can be found, in terms of the relaxation 
modulus, from the expression 

/>oc 

f]* = G{s)e~ ws ds (50) 
Jo 

Upon substituting eq|8^ for the memory function G(s) into eq ^0], one obtains the predictions of the Gaussian 
approximation for 77' and 77" . These are given by eqs |8^ in appendix |^. 

The zero shear rate viscosity r) p ^ and the zero shear rate first normal stress difference coefficient ^1,0, 
can be obtained from the complex viscosity in the limit of vanishing frequency, 

r/p.o = Urn r)'(u>) ; *i. = hm 2?1 ^ (51) 

The predictions of the zero shear rate viscometric functions by the Gaussian approximation are given by 
eqs |9(] and |9l| in appendix ||. They are compared with the exact results, eqs |J and |25|, evaluated by 
Brownian dynamics simulations, in section 7 below. 



6 First Order Perturbation Expansion 

The retarded motion expansion, eq which was obtained by carrying out a perturbation expansion of the 
distribution function ip, in terms of velocity gradients, is valid for arbitrary strength of the excluded vol 
interaction. In this section, using arguments similar to those in the papers by Ottinger and co-workers 
we derive a perturbation expansion of r p in the strength of excluded volume interaction, which is valid for 
arbitrary shear rates. A significant benefit of the perturbation expansion will be a better understanding of 
the nature of the Gaussian approximation. 

The distribution function ip may be written as ipR + ip z * , where ipR is the distribution function in the 
absence of excluded volume, i.e. in the Rouse model, and ip z * is the correction to first order in the strength 
of the excluded volume interaction. Since ipR IS Gaussian, it has the form given by eq |3^, with M{t) replaced 
by j\ffj(t), and Ojk replaced by ctX. = (QjQk/p- The second moments (QjQk) can then be expanded to 
first order as, (QjQ k ) = 0% + (Q-jQk),,- 

On substituting this expansion into eq ^, and equating terms of like order, the second moment equation 
can be separated into two equations, a zeroth-order equation and a first-order equation. The zeroth-order 
equation, which is the second moment equation of the Rouse model, is linear in c&, and has the following 
explicit solution, 



■ — (t-s)A 



7[i](M) (52) 
i'fc J 



where, £ is an exponential operator. Properties of exponential operators that operate on (N — l) 2 x (N — l) 2 
matrices are discussed in appendix [b|. The exponential operators used in this section have similar properties, 
but operate on (N — 1) x (JV — 1) matrices. 

The first-order second moment equation has the form, 

d {Q 3 Q k ) z , = « • {Q 3 Qk) z , + (Q 3 Q k ) z , ■ k t 



At 



7 E { (QjQm) z . A mk + A jm {Q m Q k ) z , } + Y% (53) 
^ ' m=l 
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where, Y R k is given by eq [l^, with the averages on the right hand side evaluated with ipR, i.e., ( ■ • • ) are 
replaced with {' ' ') R - Since ipR is a Gaussian distribution, the decomposition result, eq with ( ■ • ■ ) 
replaced with ( • • • ) , can be used to reduce Y R k to a function of second moments alone. This leads to, 



JV-l 



R 
km 



(54) 



m—1 



In eq Aj^, is given by eq |4^, with oj^ replaced by ojj. in the definition of cr^ on the right hand side. 
Equation |53f is a system of linear inhomogeneous ordinary differential equations, whose solution is, 



JV-l 

E 



;£ 



-j(t-s)A 



E(t, S )-Y«(s)-E T (t,s) 



H 
T 



(t - s)A 



(55) 



where, E is the displacement gradient tensorS 

It is immediately clear from eq [53] that the Gaussian approximation is exact to first order in the strength 
of excluded volume interaction. This follows from the fact that it could have also been derived by expanding 
eq ^3| to first order in z* . It will be seen later that this property of the Gaussian approximation, is helpful 
in elucidating its nature. 

The first order perturbation expansion for the stress tensor can be obtained by expanding Kramers 
expression, eqjll], to first order in z* . After reducing complex moments evaluated with the Rouse distribution 
function to second moments, the stress tensor can be shown to depend only on second moments through the 
relation, 

JV-l JV-l 

rP = -n p H ]T o% k -n v H £ {Q k Q k ) z , + Z R + (N - 1) n p k B T 1 (56) 



k=l 



k=l 



where, Z R is given by eq 39, with oj^ replaced by er^. in the definition of er M „ on the right hand side. 
Equations [52] and ^ may then be used to derive the following first order perturbation expansion for the 
stress tensor in arbitrary homogeneous flows, 

JV-l 



,k B T I 

i J — i 



■>£ 



2H , 
- — (t- s )A 



E{t, S ) ■ \ (k(s)+K T ( S )) S r 



H 



E T {t, s) + Z R + {N- 1) n p k B T 1 



(57) 



Note that Z R , the direct contribution to the stress tensor, is isotropic only in the limit d* — > 0. We now 
consider the special case of steady shear flow, and obtain the zero shear rate viscometric functions. 



6.1 Steady Shear Flow 



In order to obtain zero shear rate viscometric functions correct to first order in 2*, it is necessary to evaluate 
the time integrals in eqs [5^ and [37| and to evaluate the quantities Y R k and Z R in steady shear flow. The 
results of these calculations are given below, while the details are given in appendix ^|. 

The excluded volume contributions to the zero shear rate viscometric functions (correct to first order in 

to ^ of appendix ^ are, 



f ) obtained by setting 7 equal to zero in eqs 96 



(E) 
%,0 



H 



N 

E 



(d* 2 + S$) 



7/2 



C(0) c(l) 



' 9(1) 



(58) 
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A' 



(E) 
1.0 



X H Z * E 



7/2 



2,0 







(59) 



(60) 



where, the time constant \h = (£/4iT) has been introduced previously in appendix B|. and the quantities 

max(/i,y)- 1 



s\iv \ which occur in these functions and which were introduced earlier by Ottinger,ll3 are defined by, 



E c i 

j,k — min(/x,z/) 



(61) 



The first order perturbation theory predictions of the zero shear rate viscometric functions given above 
are compared with exact Brownian dynamics simulations and the Gaussian approximation in section 7. 
We first, however, examine the role of the parameters d* and z* in the present model, by considering the 
end-to-end vector at equilibrium in the limit of large N . 



6.2 The Equilibrium End-to-End Vector For Large Values of N 

The second moment of the end-to-end vector r at equilibrium is given by the expression, 

N-l 

( rr ) cq = E (QM cq 

j,k=l 



(62) 



For the Rouse model, a^ k 



(ksT/H) Sjk 1. One can show, from eqpq, that the first order correction to 



eq 



the second moments has the following form at equilibrium, 



E • rs Y • 



(63) 



where, the (N — l) 2 x (TV — l) 2 matrix Rjk, mn is defined by, Rjk, mn = Aj m Sk n + Sj m A^ n , and Yf k has the 
form, 

N 1 I N-l ^ 

(64) 



V R 
1 jk 



k B T x -» 

3 ^- E 



cq 



5/2 



Note that the function 9(fj,,m,n, v) has been introduced previously in appendix [b] (see eq |83|). It follows 
that the mean square end-to-end vector at equilibrium, correct to first order in z* , is given by, 



(r 2 ) = 

\ / CQ 



3k B T 



<=q h 



1 N 

(N-l) + -z* E 



(d* 2 + \n- v\) 



5/2 



(65) 



We now consider the limit of a large number of beads, N. In this limit, the sums in eq |6^ can be replaced 
by integrals. Introducing the following variables, 



/i v d* 

X ~ N' V ~ N' ~ ~7n 



and exploiting the symmetry in x and y, we obtain 

2X 3k B T 



oq H 



N { l + z*VN I d.v I 'd,, U 



o Jo (d 2 + x — y) 

x>y+c 



5/2 



(66) 



(67) 
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where, c is a cutoff parameter of order 1/N which accounts for the fact that fi ^ v. 

It is clear from eq [37] that the excluded volume corrections to the Rouse end-to-end vector are proportional 
to z* yN. As a result, the proper perturbation parameter to choose is z = z* yN, and not z* . This is a 
very well known result of the theory of polymer solutions,cTt3 and indicates that a perturbation expansion 
in z* becomes useless for long chains. 

The integrals in eq ^7] can be performed analytically. However, we are interested only in the form of 
eq which leads to a very valuable insight. Defining the quantity a, which is frequently used to represent 
the swelling of the polymer chain at equilibrium due to excluded volume effects, 



f 

< eq 



a 2 = f^f (68) 



we can see that in the limit of long chains, a = a (z, d). In other words, a depends asymptotically only on 
the parameters z and d, and not on the chain length N. We shall see later that this insight is very useful in 
understanding the results of Brownian dynamics simulations, and the Gaussian approximation. 



7 Equilibrium Swelling and Zero Shear Rate 
Viscometric Functions 

The prediction of equilibrium properties and zero shear rate viscometric functions, by Brownian dynamics 
simulations, the Gaussian approximation and the first order perturbation expansion, are compared in this 
section. Before doing so, it is appropriate to note that an equilibrium property, frequently defined in static 
theories of polymer solutions, namely, the swelling of the radius of gyration, a 2 , can be found from the 
following expression, 

, = i^U = Vpfl (6g) 

9 (r 2 ) r < 

\ 9 I eq 

because of the relation between the radius of gyration and the zero shear rate viscosity, eq [2^. Plots of a 2 in 
this section must, therefore, also be seen as plots of the ratio of the zero shear rate viscosity in the presence 
of excluded volume interactions to the zero shear rate viscosity in the Rouse model. 

Figures 1 to 3 are plots of a 2 , a 2 and (^i^/^fo) versus d* , respectively, at a constant value of z* = 0.5, 
for three different chain lengths, TV = 3, N — 6 and N = 12. The squares, circles and triangles are 
exact results of Brownian dynamics simulations for the narrow Gaussian potential, the dashed lines are the 
predictions of the Gaussian approximation, and the dot-dashed curves are the predictions of the first order 
perturbation theory. 

In the limit d* — > and for large values of d* , for all the values of chain length N, the Brownian dynamics 
simulations reveal that equilibrium and zero shear rate properties tend to Rouse model values. In the case 



of a and a 2 , this is expected because of the rigorous result, eq |33|. An immediate implication of this 
behavior is that, for chains of arbitrary but finite length, it is not fruitful to use a S- function potential to 
represent excluded volume interactions. On the other hand, the figures seem to suggest that a finite range 
of excluded volume interaction is required to cause an increase from Rouse model values. Both the first 
order perturbation theory and the Gaussian approximation predict a significant change from Rouse model 
values in the limit d* — ► 0. Iruthe case of a dumbbell model, we were able to rigorously understand the 
origin of these spurious results.13 The incorrect term-by-term integration of a series that was not uniformly 
convergent was found to be the source of the problem. Since first order perturbation theory is the basis 
for renormalisation group calculations, the invalidity of the (^-function potential, which is frequently used 
in these calculations, is at first sight worrisome. However, we shall see below that the use of a (5-function 
potential may be legitimate when both the limits, N — > oo and d* — > 0, are considered. 

Figures 1 to 3 show that there exists a threshold value of d* at which, the results of the Gaussian 
approximation and the first order perturbation theory, first agree with exact Brownian dynamics simulations. 



This is consistent with the first order perturbation theory predictions of the end-to-end vector, eq 65, and 



the viscometric functions, eqs d8^ and M, which reveal that, excluded volume corrections to the Rouse 
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Figure 1: Equilibrium swelling of the end-to-end vector versus the extent of excluded volume interaction d* , 
at a constant value of the strength of the interaction z*, for three different values of chain length, N. The 
non-dimensional parameters z* and d* characterize the narrow Gaussian potential, and are defined in eq 
The squares, circles and triangles are results of Brownian dynamics simulations, the dashed and dot-dashed 
lines are the approximate predictions of the Gaussian approximation, and the first order perturbation theory, 
respectively. The error bars in the Brownian dynamics simulations are smaller than the size of the symbols, 
and the continuous curves through the symbols are drawn for guiding the eye. 



model decrease with increasing values of d* . The threshold value of d* at which the approximations become 
accurate, increases as N increases. This is a consequence of the well known result, which was demonstrated in 
section 6, that excluded volume corrections scale as z* y/~N. Note, however, that the Gaussian approximation 
always becomes accurate at a smaller threshold value of d* than the first order perturbation theory. The 
Gaussian approximation, while being a non-perturbative approximation, is nevertheless, exact to first order 
in z* . Consequently, as mentioned earlier, it might be considered to consist of an infinite number of higher 
order terms, and can be expected to be more accurate than the results of the first order perturbation theory. 

All the results in Figures 1 to 3 are entirely consistent with the results obtained earlier with a dumbbell 
model for the polymer molecule. However, in the case of the dumbbell model, the dependence of the quality 
of the approximations on the chain length N, could not be examined. The results in Figures 1 to 3 seem 
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Figure 2: Swelling of the radius of gyration versus d* , at a constant value of z* , for three different values 
of N. Note that a 2 = Vpfi/Vpfi- The symbols are as indicated in the caption to Figure 1. The error bars 
in the Brownian dynamics simulations are smaller than the size of the symbols, and the continuous curves 
through the symbols are drawn for guiding the eye. 



to suggest that the Gaussian approximation has a rather limited validity, since for a given value of z* and 
d* , it gets progressively worse as the chain length N increases. This is in fact not a realistic picture — as is 
revealed below — when the data is reinterpreted in terms of a different set of coordinates. 

Figures 4 to 6 are plots of a 2 , a 2 and (^'i i o/vE , f ) versus d — {d*/\/N), respectively, at a constant value 
of z = z*\/N = 1.0, for three different chain lengths, TV = 6, N = 12 and N = 24. Before we discuss the 
figures, it is appropriate to make a few remarks about the choice of the variables in terms of which the data 
are displayed. Firstly, we choose z as the measure of the strength of excluded volume interaction because 
perturbation theory clearly reveals that excluded volume corrections scale as z* y/N. A constant value of z, 
as N increases, implies that z* must simultaneously decrease in order to keep the relative role of excluded 
volume interactions the same. Secondly, we choose the x-axis coordinate as d — (d*/\/N), because, as was 
shown in section 6, perturbation theory in the limit of long chains indicates that when the data is displayed 
in terms of d and z, all the curves should collapse on to a single curve as N — > oo. The parameter d may be 
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Figure 3: Ratio of the zero shear rate first normal stress difference coefficient in the presence of excluded 
volume interactions to the zero shear rate first normal stress difference coefficient in the Rouse model versus 
d*, at a constant value of z* , for three different values of N. The symbols are as indicated in the caption 
to Figure 1. The error bars in the Brownian dynamics simulations are smaller than the size of the symbols, 
and the continuous curves through the symbols are drawn for guiding the eye. 

considered to be the extent of excluded volume interaction, measured as a fraction of the unperturbed (i.e., 



We first discuss the results of exact Brownian dynamics simulations displayed in Figures 4 to 6. As in 
Figures 1 to 3, all the properties start at Rouse values at d — 0, go through a maximum as d increases, and 
then finally decrease once more towards Rouse values with the continued increase of d. However, as the chain 
length increases, the values seem to rise increasingly more rapidly from the Rouse values at d = 0, towards 
the maximum value. In other words, the slope at the origin seems to be getting steeper as N increases. 
Indeed, the trend of the data leads one to speculate that, in the limit N — > oo, the data might be singular 
at d = 0, and consequently legitimize, in this limit, the use of a 5-function excluded volume potential. This 
conclusion is of course only speculative, and needs to be established more rigorously. It has not been possible 
to examine more closely, with the help of Brownian dynamics simulations, the behavior at small values of d 
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Figure 4: Equilibrium swelling of the end-to-end vector versus the re-scaled extent of excluded volume 
interaction d = d*/>/~N, at a constant value of the re-scaled strength of the interaction z = z*^/N, for 
three different values of N. The squares, circles and triangles represent the results of Brownian dynamics 
simulations for N equal to 6, 12 and 24 beads, respectively. The continuous, dashed, and dot-dashed 
curves are the approximate predictions of the Gaussian approximation for N equal to 6, 12 and 24 beads, 
respectively. The error bars in the Brownian dynamics simulations arc smaller than the size of the symbols. 



for larger values of N, because of the excessive CPU time that is required. In terms of the non-dimensional 
time t* = (t/Ajj), for N = 24, a run for two non-dimensional time steps At* — 0.1 and At* = 0.08, required 
roughly 65 hours of CPU time on a SGI Origin2000 with a 195 MHz processor. 

When viewed in terms of z and d, the Gaussian approximation is revealed to be far more satisfactory 
than appeared at first sight in Figures 1 to 3. Indeed, for relatively small values of d, where the Gaus- 
sian approximation is inaccurate at small values of chain length, the Gaussian approximation seems to be 
becoming more accurate as N increases. One might expect that as N — > oo, the Gaussian approximation 
becomes accurate for an increasingly larger range of d values. However, as will perhaps become clearer with 
the results discussed below, it appears that, for a given value of z, there exists a threshold value of d, below 
which the Gaussian approximation will be inaccurate, no matter how large a choice of N is made. The 
reason for this behavior is related to a feature that is just noticeable in these figures — curves for different 
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Figure 5: Swelling of the radius of gyration versus d, at a constant value of z, for three different values 
of N. The symbols are as indicated in the caption to Figure 4. The error bars in the Brownian dynamics 
simulations are smaller than the size of the symbols. 



values of iV appear to be converging to an asymptote. This feature will become much clearer in Figure 7, 
and will be discussed in greater detail below. 

For the sake of clarity, the predictions of the first order perturbation theory are not displayed in Figures 4 
to 6. In contrast to the situation in Figures 1 to 3, where the accuracy of the first order perturbation theory 
becomes progressively worse as N increases, its accuracy appears frozen when viewed in terms of z and d. In 
other words, for different — sufficiently large — values of N, the first order perturbation theory first becomes 
accurate at the same threshold value of d. As in the case of the predictions of the Gaussian approximation, 
curves for different values of N appear to be converging to a common asymptote. This can be seen clearly 
in Figure 7. 

Figure 7 displays plots of a 2 versus d, for different chain lengths, at a constant value of z = 1. It 
clearly reveals the fact that, both in the Gaussian approximation and in the first order perturbation theory, 
curves for different values of N collapse on to a single curve in the limit N — > oo. A similar approach to an 
asymptotic limit is observed as N — > oo, in the predictions of a 2 and (^1,0/^1^0) by both the approximations, 
when they are plotted versus d. The results of Brownian dynamics simulations for N = 24 are also plotted in 
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Figure 6: Ratio of the zero shear rate first normal stress difference coefficient in the presence of excluded 
volume interactions to the zero shear rate first normal stress difference coefficient in the Rouse model versus 
d, at a constant value of z, for three different values of N. The symbols are as indicated in the caption to 
Figure 4. 



Figure 7. They indicate that for z — I , asymptotic values have already been reached by Brownian dynamics 
simulations, at this relatively small value of N, for d > 0.3. One expects that as N increases, asymptotic 
values will be reached for smaller and smaller values of d. 

The asymptotic values predicted by the first order perturbation theory were obtained by carrying out 
the integrals in eq ^ analytically. It is worth noting that the convergence to the asymptotic value is quite 
slow as d — ► 0. On the other hand, the asymptotic values predicted by the Gaussian approximation were 
obtained by a numerical procedure, as discussed below. 

In the Gaussian approximation, calculation of the equilibrium and zero shear rate quantities requires the 
evaluation of the equilibrium moments fjk- These are found here, as mentioned in appendix |b| by numerical 
integration of the system of ordinary differential equations, eq [77[ using a simple Euler scheme, until steady 
state is reached (the symmetry in j and k is used to reduce the number of equations by a factor of two). In 
addition, the evaluation of rj Pi o and ^i : o requires the inversion of the (N — I) 2 x (N — I) 2 matrix Ajk,mn 
(see eqs pu and 91). As a result, the CPU time scales as N 6 , and makes the task of generating data for 
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Figure 7: Equilibrium swelling of the end-to-end vector versus d, at a constant value of z, for different N. 
The continuous and dotted curves are the predictions of the Gaussian approximation for iV equal to 36 and 
40 beads, respectively. The filled squares are the asymptotic predictions of the Gaussian approximation, 
obtained by numerical extrapolation of finite chain data to the limit of infinite chain length. The dashed 
and dot-dashed curves are the predictions of the first order perturbation theory for N equal to 500 and 1000 
beads, respectively. The filled diamonds are the predictions of the first order perturbation theory in the long 
chain limit, obtained by carrying out the integrals in cq|(37] analytically. The circles, with error bars, are the 
results of Brownian dynamics simulations for N = 24. 

large values of N extremely computationally intensive. We have explored the predictions of chains up to 
a maximum of N = 40, since for this value of N, a single run on the SGI Origin2000 computer required 
approximately 54 hours of CPU time. The asymptotic values in Figure 7 were obtained by the following 
procedure. For z = 1, equilibrium and zero shear rate data, consisting of property values at different pairs of 
values (d, N), were first compiled by performing a large number of runs for various values of N as a function 
of d. A specific value of d was then chosen, and assuming that the various properties were functions of 1/y/W, 
the values-fbr different N were extrapolated to the limit N — > oo using a rational function extrapolation 
algorithm.!!!] The choice_pf 1/yN is motivated by the fact that the leading correction to the integrals in 
eq|7|is of order 1/VnH 
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Figure 8: Swelling of the radius of gyration versus d, at two values of z, for three different values of N. 
The squares, circles and triangles (filled for z = 5 and empty for z = 10) represent the results of Brownian 
dynamics simulations for N equal to 6, 12 and 24 beads, respectively. The continuous, dashed, and dot- 
dashed curves are the approximate predictions of the Gaussian approximation for N equal to 6, 12 and 24 
beads, respectively. The error bars in the Brownian dynamics simulations are smaller than the size of the 
symbols. 



The quality of Gaussian approximation as a function of the variable z, for the quantities a 2 and 
(^i^o/^f o)i i s displayed in Figures 8 and 9, respectively The behavior of a 2 has not been displayed as 
it is very similar to that of a 2 . It is clear from these figures that for a given value of N, the threshold value 
of d beyond which the Gaussian approximation is accurate increases as z increases. On the other hand, 
as in the case of z = 1, for a fixed value of z, the accuracy of the Gaussian approximation seems to be 
increasing with N, for small values of d. There is, however, clearly a limit to this accuracy. As ./V becomes 
large, the results of the exact Brownian dynamics simulations and the Gaussian approximation approach 
asymptotic values, and consequently, no further change can be noticed with changing N. Figures 8 and 9 
seem to indicate that at small values of d, while the asymptotic values of Brownian dynamics simulations lie 
below the asymptotic values of the Gaussian approximation for a 2 , the opposite is true for (^l^/^fo)- A 
clearer picture would be obtained if it were possible to carry out Brownian dynamics simulations with larger 
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Figure 9: Ratio of the zero shear rate first normal stress difference coefficient in the presence of excluded 
volume interactions to the zero shear rate first normal stress difference coefficient in the Rouse model versus 
d, at two values of z, for three different values of N. The symbols are as indicated in the caption to Figure 8. 
The error bars in the Brownian dynamics simulations are smaller than the size of the symbols. 



values of N. ■— , 

Typical experimental values of z lie in the range < z < 15fi9 As we have seen above, for large enough 
values of d, the Gaussian approximation remains accurate for a significant fraction of values of z in this 
range. Since corrections to the Rouse model due to excluded volume interactions decrease with increasing 
shear rate, we can anticipate that the accuracy of the Gaussian approximation will improve as the shear 
rate increases. Furthermore, since the Gaussian approximation is extremely accurate for the treatment of 
hydrodynamic interaction effects, and since hydrodynamic interaction isJikely to be the dominant effect in 
a combined theory of hydrodynamic interaction and excluded volume,E2l it is perhaps fair to say that the 
results obtained so far clearly indicate the practical usefulness of the Gaussian approximation. 
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8 Conclusions 



The influence of excluded volume interactions on the linear viscoelastic properties of a dilute polymer solution 
has been studied with the help of a narrow Gaussian excluded volume potential that acts between pairs 
of beads in a bead-spring chain model for the polymer molecule. Exact predictions of the model have 
been obtained by carrying out Brownian dynamics simulations, and approximate predictions have been 
obtained by two methods — firstly, by carrying out a first order perturbation expansion in the strength of 
excluded volume interaction, and secondly, by introducing a Gaussian approximation for the configurational 
distribution function. 

The most appropriate way to represent the results of model calculations has been shown to be in terms of 
a suitably normalized strength of excluded volume interaction z, and a suitably normalized extent of excluded 
volume interaction d. When the results are viewed in terms of these variables, the following conclusions can 
be drawn: 

1. The use of a S- function excluded volume potential (which is the narrow Gaussian excluded volume 
potential in the limit d — > 0) is not fruitful for chains with an arbitrary, but finite, number of beads 
N, because it leads to the prediction of properties identical to the Rouse model. The narrow Gaussian 
potential with a finite, non-zero, extent of interaction d, on the other hand, causes a swelling of the 
polymer chain at equilibrium, and an increase in the zero shear rate properties from their Rouse model 
values. 

2. Curves for different — but sufficiently large — values of chain length N, collapse on to a unique asymptotic 
curve in the limit N — > oo. The manner in which the results of Brownian dynamics simulations approach 
the asymptotic behavior indicates that there might be a singularity at d = 0, and consequently, the 
use of a (5-function excluded volume potential might be justified in the limit of infinite chain length. 

3. The accuracy of the first order perturbation expansion becomes independent of TV for large N. For 
a given value of z, there exists a threshold value of d beyond which the results of the first order 
perturbation theory agree with the exact results of Brownian dynamics simulations. 

4. As in the case of the first order perturbation expansion, there exists a threshold value of d beyond 
which the results of the Gaussian approximation agree with exact results. For a given value of z, this 
threshold value for accuracy is smaller than the threshold value in the first order perturbation theory. 
The accuracy of the Gaussian approximation decreases with increasing values of z. 

Explicit expressions for the end-to-end vector and the viscometric functions in terms of the model param- 
eters, obtained by carrying out a first order perturbation expansion, enable one to understand the behavior 
of the Gaussian approximation. This is because the Gaussian approximation is shown here to be exact to 
first order in z. 

The accuracy of the Gaussian approximation, for a given value of z and d, is expected to improve as the 
shear rate increases. This follows from the fact that corrections to the Rouse model, due to excluded volume 
interactions, decrease with increasing shear rate. Viscometric functions at non-zero shear rate predicted by 
Brownian dynamics simulations, the Gaussian approximation, and the first order perturbation expansion, 
will be compared in a subsequent publication. 

The advantage of using the narrow Gaussian potential to represent excluded volume interactions is 
that the accuracy of approximate solutions can be assessed by comparison with exact results. This is 
in contrast with the situation for approximate renormalisation group calculations based on a (5-function 
excluded volume potential, whose accuracy can only be judged by comparison with experiment, or with 
Monte Carlo simulations based on a different excluded volume potential. 

The results obtained here indicate that the Gaussian approximation is an accurate approximation for de- 
scribing excluded volume interactions, albeit within a certain range of parameter values. Since the usefulness 
of the Gaussian approximation has already been established for the treatment of hydrodynamic interactions, 
it is clearly worthwhile to examine the quality of the Gaussian approximation in a model for the combined 
effects of hydrodynamic interaction and excluded volume. 
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A P e q{r vix ) in the limit d — ► or d — > oo 

Upon substituting eq |lj and the Fourier representation of a <5-function into eq [}0|, and rearranging terms, 
one obtains, 



Peq( r fij,) 



■ Jdk e^-fc /jVeq Jd Ql . . . dQ N _, e~ [^ T ME::: Q,)*] 



(70) 



We now consider the integral within braces on the right hand side of eq |70j, and take up the integration over 
the bead connector vector Q 1 . Separating out all the terms containing the vector Q 1 , we can rewrite this 
integral as, 

r r h ( v ~ 1 \ 

K q dQ 2 ... dQ N _ l exp -^rj; ^Qr' D 1 " S ^ ®J k 



1 N if/" 

2k^r ^ E <yTa ~ r ^ | J dQl exp 



i,/3 = 2 
^13 



H 

2k B T 



knT 



[E{r 1 -r 2 ) + E{r 1 -r 3 ) 



Qi — iSiftQx ■ k 



■E{n-r N )] 



(71) 



where, a typical term of the excluded volume potential contribution to the Q x integral, has the form, 

v k T I 1 f 
E{rx-rp) = | 3 expj-- 5 (o; + 2g 1 T /? 2 + Oi + 2q 2 T /?3 + ... 
[2nd' ! \2 I 2ar \ 



+ 2 Qp-2 ' r 0, 0-1 + Q/3-1 



We now convert the Q x integral into spherical coordinates. In order to do so, we need to choose a reference 
vector to fix a direction in space. In the Q 1 integration, all the other vectors, Q 2 , ■ ■ ■ , Qn-x an d k are fixed. 
Without loss of generality, we choose the fixed vector as Q 2 , denote its direction as the z direction, and 
choose, in the plane perpendicular to Q 2 , an arbitrary pair of orthogonal directions as the x and y axes. Let, 
0i> 0/32 ? and 9k represent the angles that the vectors Q 1 , rp 2 and k make with the z direction, respectively. 
Similarly, let </>i, (j)^, and <f>k represent the angles that the projections of these vectors on the xy plane, make 
with the x direction. Then, 

Qx -rp 2 = Qlf/32^/32 (0l,((>l) 

where, Q\ and rg 2 represent the magnitudes of Q 1 and rp 2 , respectively, and, 



Fp2 (8i,<j)i) = sin 0i sin (cos 4>\ cos ^2 + sin^i sin ^2) + cos6>i cosi 
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Defining the function Fk(9x,4>x) similarly, we can rewrite the Q 1 integral in expression |7l|, in terms of 
spherical coordinates as, 
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dQi / dOx 

n JO JO 



•>i Q 1 sin^i exp 
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2k B T 



Ql-i5 llJt QikF k {Qx,<t>x) 
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For Qi = 0, the integrand is identically zero. For Qi ^ 0, in the limit d — > or d — > oo, the integrand tends 
to, 

— 7^— =Qf — iSifiQikF/. (01,01} 

The integrand is also a bounded function of <5i for all values d. 

An argument similar to the one above can be carried out for each of the remaining integrations over 
Q 2 , ■ ■ ■ , Qn-i- ^ follows that, 



lim / dQ l . . . dQ N _i exp 



/ H 



( 4 



dQ x . . . dQ N _ 1 exp 



\2k B T 



\k B T 
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With regard to the normalization factor A4q, since, 



dQi . . . dQ N _! exp 



we can show, by adopting a procedure similar to that above that, 



lim A4q = AC 



(73) 

(74) 
(75) 



As a result, we have established that 

1 



lim Peq(rv) 



(76) 



B Codeformational Memory-Integral Expansion 

Upon expanding the tensors ojfc in the manner displayed in eq substituting the expansion into the second 
moment equation, eq and separating the resultant equation into equations for each order in the velocity 
gradient, the following two equations are obtained, 
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with the quantities /„„ given by, 
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First Order: 



where, 
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with the quantities mn given by, 
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The function 9(n, to, n, v) has been introduced previously in the treatment of hydrodynamic interaction.!! It 
is unity if to and n lie between /x and v, and zero otherwise, 



6(jjL, rn, n, ^) 



1 if /i < to, n < or v < m,n < /i 
otherwise 



(83) 



Introducing new indices for the pairs (j,k) and (to, n), the quantity Ajk, m n may be considered an (N 
l) 2 x (AT — l) 2 matrix. The inverse can then be defined in the following manner, 



N-l 



(84) 



where, 1 is the (N — l) 2 x (N — l) 2 unit matrix IL^ „ m = <5j m (5^,„. 

In the equilibrium second moment equation, eq [f^, the term (dfjk/dt) on the left hand side, is identically 
zero. It is retained here, however, to indicate that the equation is solved for fjk by numerical integration of 
the ODE's until steady state is reached. 

Upon integrating eq |8C| with respect to time, and substituting the solution into eq |38], we finally obtain 
the expression, eq |47j, for the first order codeformational memory-integral expansion, where, the memory 
function G(t) is given by, 

N-l N-l 

G(t) = E E & m ^ & £ 

j,k=l m,n=l 

In eq 85, Xh = (C/4ff) is the familiar time constant, Tijk is defined by, 
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and the quantity Aj kj mn is given by, 
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The exponential operator £ [M] maps one matrix into another according to: 

N-l 
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S [ M ] mn ^jk. ran ~l~ Mjk, mn ~t~ T^j ^ ^ ^jk. rs ^rs, mn -\- . . . 



r,s — 1 

and has the useful properties, 
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N-l 

2J £ [ a M ] jk, rs £ [ b 1 ] rs, mn = £ [ d M + b 1] jk, mn 
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for arbitrary constants a and b. 

As mentioned in section 5, once the memory function G(s) is obtained, one can obtain the material 
functions in small amplitude oscillatory shear flow, and the zero shear rate viscomctric functions. Following 
the procedure outlined in section 5, we can show that, 

N-l N-l N-l 
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7V-1 N-l 
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Using the relations between the zero shear rate viscometric functions and rj' and rj" (eqs ^T|), one can show 
that, 

N-l N-l 

V P ,o = 4A H E E ^ikA jk rnn f mn (90) 
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N-l N-l N-l 
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C Viscometric functions correct to first order in z* 

The first step in calculating the first order excluded volume corrections to the Rouse viscometric functions, 
as mentioned earlier, is to evaluate the time integrals in eqs |5^ and ^?]. These time integrals can be carried 
out by using the forms of the tensors 7^1 and E in steady shear flow, tabulated in reference 12. One can 
show that the expression for the second moment o^j, , which is required to explicitly evaluate all the averages 
carried out with the Rouse distribution function t/jr, is given by, 

°ffc = ^jf {Sjk 1 + 2X H C jk (k + k t ) + 8X 2 H C% (k • n T ) } (92) 
while the stress tensor in steady shear flow has the form, 

N-l 

r p = -n p k B T [2A ff Cjj (k + k t ) + 8X 2 H C 2 J3 (k ■ k t )] 
j=i 
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Table 1: Functions, appearing in eq |9J, representing the indirect excluded volume contributions to the 
stress tensor in steady shear flow. The quantity SI is defined by, £1 = A/{d* 2 + sjn)) 2 . 
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(93) 



A similar expression, without the Z R term, has been derived by Ottingc in his renormalisation group 
treatment of excluded volume effects — within the framework of polymer kinetic theory — using a (5-function 
excluded volume potential. 

The next step is to explicitly evaluate the tensors Y R k and Z R , in terms of the velocity gradient k and 



the Kramers matrix Cjy, by using eq 92 for oj^.. The resultant expressions are then substituted into eq |93], 
and after a lengthy calculation, the following explicit expression for the excluded volume contribution to the 
stress tensor, correct to first order in z*, is obtained, 
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Table 2: Functions, appearing in eq representing the direct excluded volume contributions to the stress 
tensor in steady shear flow. The quantity £1 is defined by, = 4/ (d* 2 + sffl) 2 . 
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where, the functions a^f2 and /3^t, (j = 0, 1 , . . . , 6), which represent the indirect and direct contributions 
respectively, are given in Tables f and 2, and the function e M „(7) is defined by, 



e^(7) = 1 + 



(d* 2 + S$ 



( H *2 q (0)\ 0(2) _ c(l) 



(95) 



Equation |94| for the stress tensor can then be used to find the excluded volume contributions to the viscometric 
functions, correct to first order in z* , by using the definitions in eqs 
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JV 



S® (d* 2 + S$) - S$ flW + A 2 ff 7 2 (3 fig? - 2 (97) 



1 N ( 3 ) fl( 3 ) 
*f> = ^n p fe B Tz*£- = (98) 
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These expressions have been derived earlier by Ottinger,0 in an arbitrary number of space dimensions, in 
the limit d* — > 0. It is clear from eq ^8| that the second normal stress difference coefficient is zero because 
the indirect and direct excluded volume contributions cancel each other out. When d* — > 0, however, both 
the quantities a^} and fiffi are identically zero. 

References 

[1] Rouse, P.E.; J. Chem. Phys. 1953, 21, 1272. 

[2] Ottinger, H. C. J. Chem. Phys. 1989, 90, 463. 

[3] Wedgcwood, L. E. J. Non-Newtonian Fluid Mech. 1989, 31, 127. 

[4] Zylka, W. J. Chem. Phys. 1991, 94, 4628. 

[5] Prakash, J. R.; Ottinger, H. C. J. Non-Newtonian Fluid Mech. 1997, 71, 245. 

[6] Prakash, J. R. 'The Kinetic Theory of Dilute Solutions of Flexible Polymers: Hydrodynamic Interaction'; 
In Advances in the Flow and Rheology of Non- Newtonian Fluids; Siginer, D. A.; Kee, D. De; Chhabra, 
R. P.; Eds; Rheology Series; Elsevier Science: Amsterdam, 1999. 

[7] Prakash, J. R.; Ottinger, H. C. Macromolecules 1999, 32, 2028. 

[8] Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: Oxford, 1986. 

[9] des Cloizeaux, J.; Jannink, G. Polymers in Solution, Their Modelling and Structure; Oxford University 
Press: Oxford, 1990. 

[10] Schafer, L. Excluded Volume Effects in Polymer Solutions; Springer: Berlin, 1999. 

[11] Bird, R. B.; Curtiss, C. F.; Armstrong, R. O; Hassager, O. Dynamics of Polymeric Liquids. Kinetic 
Theory, 2nd edn.; Wiley-Interscience: New York, 1987; Vol. 2. 

[12] Bird, R. B.; Armstrong, R. O; Hassager, O. Dynamics of Polymeric Liquids. Fluid Mechanics, 2nd 
edn.; Wiley-Interscience: New York, 1987; Vol. 1. 

[13] Ottinger, H. C. Stochastic Processes in Polymeric Fluids; Springer: Berlin, 1996. 

[14] Schieber, J. D.; J. Rheol. 1993, 37, 1003. 

[15] Wedgewood, L. E. Rheol. Acta 1993, 32, 405. 

[16] Ottinger, H. C; Rabin, Y. J. Non-Newtonian Fluid Mech. 1989, 33, 53. 

[17] Ottinger, H. C. Phys. Rev. 1989, A40, 2664. 

[18] Zylka, W.; Ottinger H. C. Macromolecules 1991, 24, 484. 



31 



[19] Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in FORTRAN, 
2nd edn.; Cambridge University Press: Cambridge, 1992. 

[20] The scaling of linear viscoelastic properties with molecular weight, and the frequency dependence of 
oscillatory shear flow material functions, for instance, seems to be nearly entirely determined by hydro- 
dynamic interaction effects. 



32 



